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Abstract 



We have developed a Monte Carlo event generator for non-resonant diphoton (77) production 

at hadron collisions in the framework of GRQPPA, which consistently includes additional one-jet 

I^ . production. The jet-matching method developed for initial-state jet production has been extended 

' to the final state in order to regularize the final-state QED divergence in the qg ^ yy + q process. 

j^ i A QCD/QED-mixed parton shower (PS) has been developed to complete the matching. The 

PS has the capability of enforcing hard-photon radiation, and small-Q^ photon radiations that 

are not covered by the PS are supplemented by using a fragmentation function. The generated 

Qs , events can be passed to general-purpose event generators in order to perform the simulations 

down to the hadron level. Thus, we can simulate the isolation requirements that must be applied 

in experiments at the hadron level. The simulation results are in reasonable agreement with the 

predictions from RESBOS and DIPHOX. The simulated hadron-level events can be further fed 

to detector simulations in order to investigate the detailed performance of experiments. 



^ ! 1 Introduction 

Diphoton (77) production is one of the most promising channels for the discovery of the Higgs boson 

^ I having a relatively small invariant mass (< 140 GeV/c^) at the CERN LHC. However, the measure- 

00 ' ments suffer from the large irreducible diphoton background produced by non-resonant electromag- 

CO I netic (QED) interactions [I1I2]- It is necessary to understand the properties of this background, not 

only for the search for the Higgs boson but also for detailed studies after the discovery. A precise 

understanding cannot be achieved without detailed knowledge of the photon-identification capabilities 

^^ ■ of experiments. In hadron-collision experiments, photon identification is influenced by soft hadrons 

that are additionally produced by various soft phenomena, such as initial- and final-state radiations, 

hadronization, decays, underlying events and event pile-ups, because a certain isolation condition has 

to be required in order to reduce the large contamination of tt" from hadron jets. The effects of such 

soft phenomena cannot be reliably evaluated by analytical calculations or parton-level simulations. 

Hence, it is strongly desired to provide theoretical predictions in the form of Monte Carlo (MC) event 

5— i ■ generators with which we can carry out simulations down to the hadron level. 

Cd ■ ... . . ,_^ 

■ ■ The lowest-order process for non- resonant diphoton production is very simple, as shown in Fig.[TJa). 

Despite this, the next-to-leading order (NLO) correction to this process is known to be very large [21 

0]. The large correction is predominantly due to the contribution from real radiation processes, 

illustrated in Figs. [TJb) to (d). While the gluon-radiation contribution from Fig. [IJb) should not be 

very large, the contribution from quark-radiation processes shown in Figs. [Ijc) and (d) may become 

large due to the large gluon density inside protons. It is necessary to include these processes in order 

to provide a realistic simulation. However, available MC event generators consistently supporting 

these processes are limited because of a difficulty in treating the divergences associated with photon 

radiation. Although a recent release of SHERPA [S] seems to support such photon production, no 

other MC event generators are available as far as we kno\4^. 



*E-mail: shigeru.odaka@kek.jp. 

^ A forthcoming release of HERWIG-I--I- [B] may provide a consistent simulation as an NLO event generator \7\. 
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Figure 1: Typical Feynman diagrams for non- resonant diphoton production at hadron collisions: (a) 
the lowest order, (b) a gluon radiation process, (c) a quark radiation process, and (d) another quark 
radiation process. Processes (b) and (c) have initial-state QCD divergences, while (d) has a final-state 
QED divergence. 



In this report, we describe an MC event generator for non-resonant QED diphoton production in 
hadron collisions that we have developed, and discuss its predictions. The program has been devel- 
oped in the framework of the GR@PPA event generator [51 [3] and supports the generation of radiative 
processes in Figs.[ljb) to (d), which include one additional jet (light quark or gluon) in the final state, 
together with the lowest-order process in Fig.[lja). The matrix elements (ME) of these processes are 
generated using the GRACE system [lU]. Though the event generator includes radiative processes, 
it is not fully including NLO corrections because non-divergent terms in soft/virtual corrections are 
yet to be included. In any case, the radiative processes have various divergences which we need to 
regularize. The initial-state strong-interaction (QCD) divergences can be regularized using a method 
that we have developed for weak-boson production processes, where divergent components are numer- 
ically subtracted from the matrix elements of radiative processes (the LLL subtraction) [TTJ [T^ [T3] . 
The subtracted components are restored by combining non-radiative processes to which a parton 
shower (PS) is applied. We can avoid the double-count problem by the subtraction and regularize the 
divergences as a result of the multiple radiation in PS. 

The quark-radiation processes illustrated in Figs. [Ijc) and (d) has not only an initial-state QCD 
divergence but also a QED divergence in the final state. We have extended the method for initial- 
state QCD divergences to this final-state QED divergence. The extension of the subtraction method 
is straightforward, while the preparation of an appropriate PS is not trivial since it has to support 
QED together with QCD. This PS has to be applied to the qg — )■ ^q process to radiate photons 
from the final-state quark based on the collinear approximation (the fragmentation process). We 
have successfully developed such a PS on the basis of the QCD PS included in the GR@PPA 2.8 
distributioro. 

Since parton showers are based on perturbative calculations, it is necessary to apply a certain 
cutoff to avoid collinear (small-Q^) divergences. Further soft phenomena in QCD including non- 
perturbative effects are simulated by adopting appropriate models in general-purpose event generators, 
e.^., PYTHIA [TT. However, such simulations are not available for QED photon radiation. Since 
we consider radiations from final-state partons, energetic photons visible in detectors can emerge 
even from very small-Q^ branches. In order to simulate such small-Q^ photon radiations, we have 
adopted a method employing a fragmentation function (FF) . The PS supplemented with an FF-based 
simulation that we have developed has the capability of enforcing hard-photon radiation. This function 
dramatically improves the generation efficiency for the fragmentation process. 



^The base release is the GROPPA 2.8.3 update [9]. 



The generated events can be passed to PYTHIA [14] for simulating small-Q^ QCD phenomena 
down to the hadron level. Hence, we can evaluate the effect of isolation requirements at the hadron 
level using the simulated events. Those events processed by PYTHIA can be fed to detector simulations 
of experiments. Further detailed studies will become possible with such simulations. 

We sometimes find reports in which the detection efficiency and acceptance for diphoton mea- 
surements are evaluated by using event generators for the lowest-order process and the fragmentation 
process. Such evaluations are not self-consistent. Parton showers used in the simulation of the frag- 
mentation process have an energy scale to be determined arbitrarily. The energy scale defines the 
maximum hardness of the radiation. The simulation results depend on this energy scale since those 
radiations exceeding the scale are ignored. Non-collinear contributions are also ignored. These ignored 
contributions may be small in many processes compared to the contributions taken into consideration. 
However, as we will show in this report, they may become comparable to the lowest-order contribution 
in diphoton production. It is necessary to include a simulation based on exact matrix elements for 
radiative processes in order to make a reliable evaluation. 

The fragmentation process that we take into consideration in this report is the so-called "one- 
fragmentation" process. The "two-fragmentation" process in which two photons are radiated from 
final-state partons, for instance, from quarks in gg — > qq, is not supported. Besides, the gluon-fusion 
process, gg — > 77, and its higher orders are not included in the present study. 

We require a typical kinematical condition for the Higgs-boson search at the LHC through the 
present study because we are interested in its background; that is 

Pt{ii) > 40 GeV/c, pt(72) > 25 GeV/c, 
|ry(7)| < 2.5, Ai?(77) > 0.4 
80 < rrij^ < 140 GeV/c^. (1) 

We apply an asymmetric requirement to the transverse momenta (pr) of photons with respect to 
the incident beam direction. The requirement on the pseudorapidity (ry) is common to the two 
photons. In addition, though this is not effective for real diphoton events now we consider, we require 
a sufficient AR separation between the two photons, where AR is defined from the differences in the 
pseudorapidity (77) and azimuthal angle {(f>) as AR^ = Arf + AcjP . Finally, the invariant mass of the 
two photons (fTiy-y) is restricted to the range that we are interested in. The simulations are carried 
out for the design condition of the LHC, proton-proton collisions at a center-of-mass (cm) energy of 
14 TeV, with CTEQ6L1 [TS] for the parton distribution function (PDF). 

This report is organized as follows: the extension of the limited leading-log (LLL) subtraction 
method to the final-state QED divergence is described in Section 2, and the QCD/QED-mixcd PS 
that we have developed is described in Section 3. The matching in the final-state QED radiation is 
discussed in Section 4. Results from a full simulation down to the hadron level are presented and the 
effect of a typical isolation requirement is discussed in Section 5. Physical results obtained from the 
full simulation are presented and discussed in Section 6, and the discussions are concluded in Section 
7. 



2 Final-state QED LLL subtraction 

We approximate the final-state QED divergence in the matrix element for the qg — >■ 77 -t- q process as 

■^[™?(s,*77.)f = |x,9->^,(so,$^,)f 4^^,f H0^2) 9{s,so) ^(mLl-Q')- (2) 

The leading-log (LL) radiation function can be given as 

4!;^,fHo^z) = :^^p,^,.(z) (3) 

with the splitting function that is defined as 

A^g7W = e^V-T' (4) 



27r 


Q^ 


-L q^q-i 


-4 


1 + 

1- 


Z 





where e^ is equal to 1/3 for down-type quarks and 2/3 for up-type quarks. The electromagnetic 
coupling, a, is assumed to be constant and equal to 1/132.51, the value used for matrix-element 
calculations [TB]. We evaluate Eq. ([2]) for two possible combinations of (7-7 pairs in the final state, 
and numerically subtract them from the matrix element of qg — >■ 779, together with the LLL term for 
the initial-state QCD divergence 0. The parameters Q^ and z are defined by using the sum of the 
energy {Eq^) and momentum (pq-y) of the considered q--f pair. They are defined in the cm frame of 
the qg — > 77g event as 

Q^ = EI-pI^, (5) 

and 

PL{Eq^+Pqj)+Q^/2 
PqiiEq-i +Pq-,)+Q'^ 

where pl is the momentum component of the q in parallel to the summed momentum. Equation (|6]) is 
defined so that z should represent the momentum fraction at the infinite-momentum limit, Q^ /pi^ — >■ 
0. 

Equation ^ is slightly different from the radiation function for the initial state [TT]; z is absent 
in the denominator. This factor is included in the initial-state radiation function in order to take into 
account the change of the flux factor in the cross-section calculation, which is equal to the squared cm 
energy of the considered process. The relation between the squared cm energy of the non-radiative 
subprocess (sq) and that of the radiative process (s) can be expressed as sq = zs for initial-state 
radiations, whereas final-state radiations do not alter the cm energy at the limit Q^ /p^ — > where 
theories are usually discussed. Therefore, if we naively follow this argument, we should set sq = s and 
5(s,so) = linEq. ©. 

In our PS model, the kinematics of PS branches are so defined that the momentum is strictly 
conserved whereas the energy conservation is ignored [9,. The energy conservation is restored by 
increasing the initial-state parton momenta after completing the PS simulation. The mapping from 
a gg — >■ 77(7 event to a non-radiative qg — >■ 79 event, which is necessary for evaluating the LLL term 
in Eq. ([2]), has to strictly reverse this procedure in order to achieve a good matching. We define the 
quark momentum of the non-radiative process as the sum of the momenta of the quark and one of 
the two photons in the radiative process. Hence, for non-zero Q^, sq is necessarily smaller than s. 
Accordingly, the flux factor also changes in the final-state radiation. In order to take this change into 
account, we define the correction factor in Eq. ^ as 

5(s, So) = s/so- (7) 

We note that this kind of correction, which is effective only at large Q^ , strongly depends on the 
applied PS model. The correction in the above must be appropriate only for the PS model that we 
have adopted. 

An event generator implementing the above subtraction has been developed in the framework of 
GR@PPA, and the generation was tested for the 14-TeV LHC condition. The energy scales were 
defined as [8] 

A*' = |pr(7i)-PT(72)lV4, (8) 

where pt(7i) denotes the transverse momentum vector of the two photons. This definition is equivalent 
to the ordinary definition, /i = ^7^(7), for qq — > 77 events. We used the identical definition for the 
renormalization and factorization scales. The energy scale for the initial-state PS must be equal to 
the factorization scale in our method. Though it is not necessary, we adopted the same definition for 
the final-state PS. The energy scale to limit the LLL subtraction, /illl in Eq. ([2]), must be equal to 
the energy scale of the final-state PS to be applied to the mapped non-radiative event (9j. We defined 
it to be equal to the pr of the mapped qg — ?► 79 event in the current study. 

A result is shown in Fig. [5J We have plotted the distribution of Ai? between the photon and quark 
in the final state of the 95 — )> 77 + q events, where q represents any quark or anti-quark up to the 
h quark. We obtain two values since there are two photons in the final state. We take the smaller 
one as A_R(7-jet). The constraint on the photons in Eq. ([1]) was required at the event generation. 
Parton showers are yet to be applied but the initial-state QCD LLL subtraction has been applied in 
this simulation. Therefore, the result is finite even though we have required practically nothing for 
the final-state quark. The quark is allowed to be very soft and out of the detection. This is why we 
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Figure 2: Ai?(7-jet) distribution of simulated qg — >■ 77 + 9 events satisfying the kinematical condition 
in Eq. ([IJ. The dotted histogram shows the distribution before applying the final-state QED LLL 
subtraction, where Ai?(7-jet) > 0.1 is required to cut off the divergence. The open circles and triangles 
show the distribution of positive- and negative-weight events, respectively, after the subtraction. The 
final result after applying the final-state subtraction is shown with filled circles. The initial-state 
QCD LLL subtraction has already been applied while parton showers are yet to be applied in this 
simulation. 



denote the process as qg -^ 77 + q instead of qg -^ 779- The distribution without the application of 
the final-state LLL subtraction is shown by the dotted histogram for comparison, in which we have 
applied a cutoff of Ai?(7-jet) > 0.1 in order to avoid the divergence at Ai?(7-jet) — 0. 

Since the LLL subtraction is unphysical, we obtain negative-weight events as well as ordinary 
positive- weight events when we apply the subtraction. The event weights are always equal to +1 or 
— 1 because BASES/SPRING [T71IIH] automatically unweights the events. Therefore, we can obtain 
the desired distribution by subtracting the number of negative-weight events from that of positive- 
weight events in each histogram bin. The open circles and triangles in the figure show the distribution 
of positive- and negative- weight events, respectively, and the final distribution is shown by filled circles. 
We can see that the distributions after the subtraction converge to finite values as Ai?(7-jet) — )■ 0; 
not only does the final result converge to zero, but the distributions of positive- and negative-weight 
events individually converge to a finite value. These facts imply that the subtraction was performed 
properly, at least near the divergence limit. We have applied a very small cutoff, Ai?(7-jet) > 0.01, 
in the event generation for numerical stability. We can see that the effect of this cutoff is negligible. 

The LLL subtraction is limited by the 9 function in Eq. ©. A steep rise in the final distribution 
between Ai?(7-jet) = 1.5 and 2.0 corresponds to this limitation. The subtraction is not effective in the 
large Ai?(7-jet) region where events composed of two photons widely separated from the final-state 
quark exhibit. The distribution below the rise shows a non-logarithmic contribution remaining after 
the subtraction. The distribution of negative-weight events steeply decreases at large Ai?(7-jet) since 
they emerge only in the region where the subtraction is active. The negative-weight events at very 
large Ai?(7-jet) were produced by the initial-state QCD LLL subtraction. 



3 Fragmentation process 

The subtracted LL components have to be restored by non-radiative processes to which an appropriate 
PS is apphed. The initial-state QCD LL components in gg — > 77 + 5 and qg ~> ■yy + q can be restored 
with qq — )> 77 events to which an ordinary initial-state QCD PS is applied, as in the case of weak- 
boson production. On the other hand, the final-state QED LL components in qg — ^ 77 -I- q have to be 
restored with qg — >■ jq events to which a final-state PS including QED photon radiations is applied. 
We have developed such a PS. 

3.1 QCD/QED-mixed PS 

The implementation of parton showers is very simple. First of all, we define the Sudakov form factor 
which provides the probability that no branch happens between certain energy scales. It can be 
expressed by using the radiation probability, Ti{Q'^), as 



SiQlQl)=e^p 



dQ'Y.^,{Q') 



(9) 



where the sum is taken over all relevant branch modes. The radiation probability can be obtained by 
integrating the splitting probability, ^((3^,z), over the "momentum fraction" z as 

r»(Q')= 1^^^^^ ^^{Q^,z)dz. (10) 

Although the integration should be done from to 1 ideally, we must frequently limit the range in 
order to regularize the divergences in the splitting probabilities. In our QCD PS [S], the splitting 
probability is approximated as 

1 a^'^ 
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by using the 1-loop QCD coupling, as{Q )j and the leading-order splitting functions, Pi{z). 

In the actual implementation of PS, we first determine the Q^ of a branch. Suppose that we have a 
branch at Q^ — Qp^c- In the final-state (timelike) PS, the Q^ of the next branch can be determined by 
solving the equation, S{Q^, Qp^c) — V^ by introducing a random number rj which uniformly distributes 
in the range from to 1. Once the Q^ is determined, we can determine the branch mode i according 
to the radiation probability Ti{Q'^), and then determine the z parameter in proportion to 'i>i{Q^,z). 
The thus determined PS parameters, Q^ and z, are converted to the momenta of participating partons 
based on a certain kinematics model. The model that we have adopted is discussed elsewhere [51 [HI [T^ . 

We have extended the above discussion to QED. The Sudakov form factor and the radiation 
probability for QED radiations of quarks can be defined in the same way as in the above by using the 
QED splitting probability defined as 

*QED(g',^) - ^ ^Pg^,^{z). (12) 

The splitting function Pq^q^{z) is defined in Eq. (U). We can consider this QED branch as one of the 
branch modes for defining the total Sudakov form factor in Eq. ^ , and apply the same procedure for 
determining the PS parameters, Q^ and z. Although it is not straightforward to solve the equation 
'S'(Q^, Qpro) = V because the Q^ dependence is different in the QCD and QED Sudakov form factors, 
we can obtain an accurate solution by performing an iteration. We found that if an appropriate 
iteration is applied we can achieve an accuracy better than 10^^" for the S'(Q^, Qpre) value with three 
or four iteration steps. 

The above is a straightforward extension. The so-called "old model" of the PYTHIA PS [14], 
PYSHOW, employs a different method for generating QED branches. The Q^ of the next branch is 
determined independently for possible QCD and QED branches, and the branch giving a larger Q^ 
is taken as the next one. Although the implementation is very different, these two methods give an 
identical probability for QED radiation. When we have a branch at Q^ — Qp^^, the probability that a 

6 



QED branch at Q^ — Qqed is the next branch is S{QQ^Y)TQprc)^QEJ:>{QQED) ^'^^ the both methods, 
where S'(Qqeq, Qpi,,) is the total Sudakov form factor between the two energy scales and rQED(QQED) 
is the QED radiation probability at Q^ = Qqed- We have confirmed that these two methods give 
an identical result using our simulation program. Since the result is identical, it is better to use the 
PYSHOW method because it is faster; the iteration in the first method takes time. 

A mixed PS which equally treats QCD and QED branches has been constructed by adding the 
above procedure to our PS. Incidentally, the implementation of PS is terminated at a certain Q^ 
(= Qq). We set Qo = 5 GeV in our PS. As has been discussed in the Introduction, hard photons 
can be radiated with even smaller Q^. In order to simulate such photon radiations, we add QED 
branches to partons remaining after the PS implementation according to the FF in Ref. [12]. We 
use the source code included in the DIPHOX 1.2 distribution 3 . Although this FF is based on 
the next-to-leading-log (NLL) approximation, and thus the gluons may radiate photons, we add the 
branches to the remaining quarks only. This FF program has two options. The difference between 
the two options predominantly results in different radiation probabilities from gluons. Since we do 
not consider the radiation from gluons, the difference is not significant in our application. We adopt 
the first option in the current study. 

The FF gives us the probability of the photon radiation having a certain momentum fraction of x. 
The probability is integrated over the Q^ range from a given energy scale (Qff) down to, in principle, 
Q^ = 0. We set Q|p = Qq. Though the FF gives us longitudinal momentum information, it does 
not provide any information on the transverse motion. However, this is not a serious problem in our 
application because the relevant Q^ values are always small. Although it may be sufficient to add 
9 ~^ 97 branches without any transverse motion, we add a non-zero pt to the branch products with 
respect to the parent-parton direction. The pt value is randomly chosen according to a Gaussian 
distribution. We set the standard deviation of the Gaussian distribution to 1 GeV/c tentatively in 
the present study. 

3.2 Forced photon-radiation PS 

A mixed PS supplemented with an FF radiation has been successfully constructed as described in the 
previous subsection. However, since the QED coupling is markedly smaller than the QCD coupling, 
it is not efficient to generate hard photons that we are interested in by using such a PS. It is desired 
to introduce a mechanism to enforce hard-photon radiation. 

First of all, we consider the probability that the first QED branch is at Q^ — Qqed- The probability 
of there being a QCD branch at Q^ — Qp^o is rQcoCQprc)- The probability of a QED branch 
at Q^ = Qqed being the next branch is ^'(Qqed, Qp,.o)rQED(<5QED)' ^s described in the previous 
subsection. Of course, there should be no QED branch from a given PS energy scale (/i^) down to 
Qpie- Therefore, the probability that we want to evaluate can be expressed as 

2 

+ / ^ rfQpro'5'QED(Qprc7A*^) X rQCD(Qpro) ^ 'S'(QqeD J '9prc)rQED (Qqed) • (13) 

We have added the first term in order to include the probability that the QED branch is the very first 
branch. Since S{QlQl) ^ Sqcd{QIQI)Sqeb{QIQI), SiQlQ^)S{Q^Ql) ^ S{QlQl) and' 

4^s,iQl Ql) - -sm\.Ql)'^^{Ql). (14) 

the expression in Eq. (J13p can be converted to 

^QEd(QqED'M^) = 'S'QED((5|ED'A^^)rQED(QQED)- (15) 

This is the probability of finding the first QED branch at Q^ — Qqed without taking QCD branches 
into consideration. Namely, the Q^ of the first QED branch can be determined independently of QCD. 
We have verified this conclusion using our simulation. 

7 
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A PS which enforces hard-photon radiation has been constructed based on the above discussion. 
We consider only one photon-radiation from each prompt quark produced by the hard interaction, and 
do not take into account the radiation from secondary quarks which are produced from gluon sphtting. 
This assumption must be sufficient because the photon-radiation probabihty is very small. The FF 
photon radiation is also considered for the primary quarks only. In addition, we apply a constraint 
that the radiated photons always have large energy (> -Emin)- The requirement of a hard-photon 
radiation results in certain constraints on random numbers used for executing the PS. The constraints 
determine the weight of the generated event. If this event weight is fed to the BASES/SPRING 
system appropriately, BASES/SPRING can automatically perform the cross-section integration and 
event generation according to the constraints that we require. The actual implementation of our PS 
is described in the following. The application of this PS improved the event-generation efficiency by 
a factor of 88 in the study to be described later. 

Before starting the event generation, we evaluate the total FF-radiation probability for each quark 
flavor. We can evaluate it since Qff is always equal to the Qq of PS. When a quark to which we should 
apply the PS is selected, we can evaluate the total PS photon-radiation probability from the Sudakov 
form factor Sqei){Qo, M^)- The selection of the quark is trivial in the present study because each event 
includes only one quark in the final state. We estimate the total photon-radiation probability by adding 
the FF and PS probabilities, and take this probability as the base event-weight. Then, we determine 
whether the photon should be radiated from FF or PS according to the evaluated probabilities. If 
PS is chosen, we determine Qqed at which the QED branch should happen by solving the equation 
'S'qed(Qqedi A*^) ~ V- The random number t] is constrained in the range to give a solution within 
Qq < Qqed < M^i because this constraint is already taken into account in the event weight. 

The QCD PS is applied as usual after the above procedure. During the evolution, when Q^ for 
the next branch becomes smaller than Qqedj ^ QED branch is inserted instead of a QCD branch. Of 
course, the insertion is done only if the PS radiation is selected. Then, we determine the z parameter 
of the QED branch. We apply a hard-photon constraint here. Though the full three-dimensional 
kinematics is yet to be determined, the energy of the branching parton is available at this stage of 
the simulation. Thus, the requirement of E-y > E^min can be converted to an allowed range of the z 
parameter. This constraint results in an additional event weight. The weight is zero if the energy of 
the branching parton is not large enough. After completing the determination of the z parameter, the 
QCD evolution is restarted by setting the maximum Q^ to Qqed- After completing the PS simulation, 
a photon radiation is added according to the FF if the FF radiation is selected. The hard-photon 
requirement gives another event weight here because it constrains the allowed range of the x parameter 
in the FF. 

The total event weight is evaluated as the product of all event weights, and passed to BASES/SPRING 
using the LabCut framework which is experimentally implemented in GR@PPA 2.8. Parton show- 
ers can be applied in both the integration and event-generation stages of BASES/SPRING in this 
framework. A cut can be applied to the events after the PS simulations. This framework is useful 
for generating events in which the kinematics of the particles of interest are strongly affected by PS, 
such as photon production with relatively small-p^, because in principle we can apply a very loose 
condition at the hard-interaction generation. Incidentally, no cut is applied in the present study; the 
LabCut framework is used only for passing the event weight to BASES/SPRING. We have confirmed 
that the constructed PS (the forced PS) gives results identical to those from the mixed PS above the 
hard-photon requirement, E^ > -Emin- 

3.3 Corrections 

The kinematics model which determines the relation between the PS parameters, Q^ and 2, and the 
parton momenta in the PS branch is described elsewhere [Hlin]- The model is designed so that the 
application of the final-state PS should never alter the momenta of the other particles in the final 
state. The application of a final-state PS necessarily increases the total energy in the final state since 
it makes on-shell partons develop into massive jets. The energy conservation is restored by increasing 
the initial-state parton momenta in our model. In the course of the present study, we found that 
this model is too naive to properly simulate the particle momenta inside jets. The implementation 
of PS is theoretically justified only at the collinear limit where the increase of the energy can be 



ignored. The implementation may substantially depend on the applied model away from the limit. 
Our PS model increases the squared cm energy of the hard process, s, and the cross sections of the 
hard processes that we are now interested in commonly decrease in proportion to 1/s. Hence, it must 
be natural to suppress the PS implementation in proportion to 1/s. In addition, since we increase 
the initial-state parton momenta, we should also take the change in the parton distribution function 
(PDF) into consideration. 

In order to take these s effects into account, we apply a correction to the generated events as 
follows. The correction factors are evaluate as event weights. The weight relevant to the cross-section 
change is defined as 

Ws = s/s', (16) 

where s and s' represent the squared cm energy before and after the application of the final-state 
PS, respectively. The initial-state parton momenta are increased so that this increase should not 
alter the rapidity of the hard-interaction system. Hence, the momentum fraction of the initial-state 
partons is adjusted as x'^ — ^J s' / sxi where i represents two initial-state partons (1 or 2). The weight 
corresponding to the change of the parton density is defined as 

WPDF = -pT 2TT~i 2T' \^'> 

where fi{xi,^%) is the PDF for the parton i. We assume that the factorization scale, ^%, is unchanged. 
The total weight is evaluated as w = wgWpDF. Since this weight is always smaller than one, we can 
apply a correction with a simple hit/miss rejection. If the event is rejected, the final-state PS is retried 
until the event is accepted. Therefore, this correction never alters the cross section of the event nor 
the momenta of final-state particles to which the PS is not applied. The correction solely alters the 
momentum distribution of particles inside PS and the resultant jet-mass distribution. Of course, the 
correction vanishes at the collinear limit. 

3.4 Simulation result 

The Ai?(7-jet) distribution of the fragmentation process simulated with the developed PS is shown 
in Fig. [31 The forced photon-radiation PS described in the above has been applied to the qg — )■ jq 
events generated in the LHC condition. The PS energy scale, (J-fsr, was chosen to be equal to the px 
of the generated events. The generation condition was sufficiently relaxed in order not to affect the 
final result. The initial-state PS is yet to be applied. We required a constraint of E^ > 20 GeV in 
the forced PS. The kinematical condition in Eq. ([IJ was required after the PS simulation. Thus, the 
constraint in the forced PS does not affect the obtained result. 

In order to derive the quantity Ai?(7-jet), we need to determine the momentum of the remnant 
jet in which the detected fragmentation photon is excluded. We define it from the initial-state parton 
momenta and the observed two photon momenta, without using the momenta of partons produced 
by PS. Although the obtained jet momentum is equal to the summed momentum of remnant partons 
in our simulation, this definition can be universally applied to other simulations in which the origin 
of final-state partons may be ambiguous. In any case, by using the obtained remnant-jet momentum, 
we can reconstruct a qg -^ ^^ + q event which can be compared with the parton-level qg ^ ^^ + q 
events generated according to the matrix elements. Since the jet mass of the remnant jet is ignored, 
the cm energy of the reconstructed qg ^ "J^y + q events is smaller than that of the simulated events. 
The energy conservation is restored by adjusting the initial-state parton momenta as is done in the 
PS application. 

In Fig. [21 together with the total yield (plots), we have illustrated the PS (solid histogram) and 
FF (dashed histogram) contributions separately. The FF contribution is about 40% of the total, and 
is concentrated at small Ai?(7-jet). An unnatural bumpy structure at small Ai?(7-jet) evident in 
the total distribution suggests that the transverse motion added to FF branches must be too small. 
However, we do not need to be greatly concerned about this detail because Ai?(7-jet) is a quantity 
which is hard to accurately measure in actual experiments, especially at small A_R(7-jet) where the 
fragmentation photon is confined inside a hadron jet. The tail in the FF contribution can be attributed 
to an effect of multiple QCD branches in PS. The distribution at large Ai?(7-jet) is dominated by 
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Figure 3: Ai?,(7-jet) distribution of the fragmentation process. The forced PS with a constraint of 
E^ > 20 GeV has been apphed to the qg — >■ 75 events generated in the LHC condition. The initial-state 
PS is not apphed. The kinematical condition in Eq. ([Ij is required after the PS simulation. Together 
with the total yield (plots), the PS (solid histogram) and FF (dashed histogram) contributions are 
separately shown. 



the PS contribution, 
limitation by /ipsR- 



The distribution is strongly suppressed above Ai?(7-jet) = 1.5 owing to the 



4 Matching test 

The simulation of the qg — > 77 + q process is completed by adding the two simulations, the LLL- 
subtracted qg —> 'J'J + q process and the fragmentation process. The matching between the two 
simulations can be rigorously tested at the parton level by using the distribution of a parameter 
Rq = Q/pto- In the fragmentation events, Q is defined as the square root of the Q^ parameter 
of the PS branch in which the fragmentation photon has been produced, and the px of the prompt 
photon is taken as pro- Thus, Q is not an observable quantity. Observable kinematical variables are 
affected by multiple QCD radiations in the fragmentation events. Such an alternation is a higher-order 
effect which is not taken into account in the parton-level 95 — >■ 77 + <J events. We can perform an 
unambiguous test at a common perturbation order by using the thus-defined quantity Q. 

In our PS model, the parameter Q is assumed to be the invariant mass of the branch products. 
Hence, in the LLL-subtracted qg —> Jj + q events, the quantity that should be compared with Q is the 
invariant mass of a j-q pair. There are two possible values since two photons exist in the final state. 
We take the smaller one as Q, and take the pt of the unpaired photon as pro- This pt is the same 
as Pt of the non-radiative subprocess assumed in the LLL subtraction. Because we take the px of 
the non-radiative process as the final-state PS energy scale, fiFSR, in the fragmentation events, and as 
the upper limit of the subtraction, /xlll, in the LLL-subtracted events by default, the LL component 
must be sharply shared between the two simulations at Rq = 1 in the default setting. 

The simulation results are shown in Fig. SI The distributions were obtained with the simulations 
described in previous sections. The distribution of the fragmentation events is shown with a solid 
histogram, and that of the LLL-subtracted qg — ^ 77 -I- q events with a dashed histogram. The dashed 
histogram at Rq < 1 shows the remaining non-LL component in the qg -^ ^^ + q process. We set 
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Figure 4: i?Q distributions of simulated fragmentation events (solid histogram) and LLL-subtracted 
ggf — >■ 77 + g events (dashed histogram). The sum of the two distributions is illustrated as a plot. 
See the text for the definition of the variable. A naive prediction from the qg — > 77 + q ME with a 
constraint of Ai?(7— jet) > 0.1 is shown with a dotted histogram for comparison. 

Q = for FF branches. Hence, the FF-radiation events in the fragmentation process are concentrated 
in the first bin. The distribution for the fragmentation process that is extended up to Rq = 1 is 
determined by PS branches. 

The two simulation results have sharp edges at Rq = 1 as expected. The sum of the two distribu- 
tions is plotted with filled circles. We can see that the two simulations are connected smoothly at the 
boundary, and that the summed distribution has a natural shape without any remarkable discontinu- 
ity. The naive prediction from the qg — >■ 77 -I- q ME with a constraint of Ai?(7-jet) > 0.1 is drawn 
with a dotted histogram for comparison. The summed distribution is in good agreement with this 
prediction down to Rq ^ 0.3. The simulation result has to be suppressed with respect to the naive 
ME-prediction at small Rq since the former is finite whereas the latter is divergent. We find that the 
suppression is effective only below Rq ^ 0.3. 

The agreement between the simulation and the ME prediction at Rq > 1 is trivial since no 
subtraction is applied in this region. On the other hand, the agreement at Rq < 1 is not trivial 
because the fragmentation result, which dominates the distribution in this region, depends on the 
applied corrections. The contribution from the non-LL component is also significant. The marked 
difference between the PS distribution and the ME prediction is completely supplemented by the 
non-LL contribution. We have confirmed that we always find a significant mismatch if any of the 
corrections described in previous sections is disregarded. Incidentally, the matching is not perfect. A 
small enhancement with respect to the ME prediction can be seen in the Rq range from 0.3 to 0.7. 
The enhancement is at most 10% and far smaller than the higher-order QCD effects to be described 
later. 

We set /ipsR = fJ-LLL = Pto in the above study. This choice is arbitrary. If the matching is perfect, 
the result should be unchanged even if we take other settings. We have tested the stability by changing 
the definition as /i = /ipsR = Mlll = O.SpTo and 1.5pTo- The boundary in the Rq distribution moves 
according to this change, as shown in Fig. [5] The sum of the distributions is over-plotted in Fig. IH] 
for the three settings. Despite the fact that the small enhancement at medium Rq that we found 
in Fig. 21 produces a small mismatch in the case oi fi = O.bpTo, the three distributions are hard to 
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Figure 5: Energy-scale dependence of the Rq distributions for the fragmentation process (a) and the 
LLL-subtracted 95 — > 77 + g process (b). The resuhs are shown for three settings: fi/pro — 0.5 
(dashed), 1.0 (sohd) and 1.5 (dot-dashed), where /i — /ipsR = Mlll. The naive ME-prediction with 
Ai?(7— jet) > 0.1 is shown with a dotted histogram for comparison. 
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Figure 6: Energy-scale dependence of the Rq distribution. The sum of the distributions in Fig. [S] is 
plotted for three settings for the energy scale: n/pro = 0.5 (open circles), 1.0 (filled circles) and 1.5 
(open triangles). The dotted histogram shows the naive ME-prediction with Ai?(7— jet) > 0.1. 
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Figure 7: As for Fig. |H but Q for the fragmentation events is defined as the invariant mass between 
the photon and quark (m^j) in the reconstructed parton-level qg ^ JJ + q events. 

distinguish from each other. Namely, the distribution is very stable against the change of the energy 
scale. The summed cross section is also stable within ±6%. 

The matching study described above has almost nothing to do with the kinematics model in 
the PS. The model only indirectly contributes to the tested distribution through the s corrections. 
Figure [7] shows another Rq distribution, in which Rq for the fragmentation events is defined using the 
reconstructed parton-level qg -^ ^^ + q events in the same way as in the LLL-subtracted qg — )■ 77 + q 
events. Hence, the distribution should reflect the effect of the kinematics model. As we can see in 
Fig. El the change of the definition does not very significantly alter the distribution. The distribution 
for the LLL subtracted gg — > 77 + g is the same as that shown in Fig.^l The FF radiations contribute 
to the distribution at small Rq in this result. 

The matching is still good, although we can see a significant enhancement in a wide range of Rq 
with respect to the naive ME-prediction, and a sight discontinuity at Rq = 1. These alternations can 
be attributed to the effects of QCD radiations preceding the QED radiation. The previously used Q 
must be identical to the 75 invariant mass in our model if the QED branch is the first branch in PS. 
Preceding QCD branches smears this identity. The smearing results in the enhancement in the Rq 
distribution because the distribution is a decreasing function of Rq. Although they are not visible 
in Fig. [3 a very small but substantial amount of fragmentation events exhibit even at Rq > 1 as 
a result of the smearing. This leakage causes the discontinuity at Rq — 1. Since the smearing is 
a higher-order effect which is not included in the lowest-order qg ^ ^^ + q ME, these alternations 
should be considered as an advantage of our simulation. A similar enhancement in the weak-boson 
Pt spectrum at medium pT is important to precisely reproduce the measurement data 13 . Note that 
the final-state QCD radiation effects are not taken into account in RESBOS and DIPHOX. 



5 Full simulation and isolation cut 

It is necessary in actual experiments to require an isolation condition in the identification of photons in 
order to reduce the huge background from hadron jets. It is difficult to reproduce the cuts applied by 
experiments with parton-level simulations. This is the main reason why hadron-level event generators 
are desired. Besides, if hadron-level events are generated consistently, the generated events can be 
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Figure 8: Ai?(7-jet) distribution of the fully simulated radiative processes. The fragmentation process 
(solid histogram) and the LLL-subtracted (?<?—>■ 77 + g process (dashed histogram) are combined to 
a smooth distribution (fulled circles), which should be compared with the naive qg — >■ 77 + q ME- 
prediction (dotted histogram). The dot-dashed histogram shows the distribution of qq — > 77 + g 
events. 

passed to detector simulations for more detailed studies. 

In this section, we simulate a typical isolation condition by using events which have been simulated 
down to the hadron level. The event generation was carried out in the 14-TeV LHC condition sepa- 
rately for the fragmentation process and the other processes including the qq — >■ 77, qq ^ JJ + g, and 
qg -> 77 -I- g processes, where q denotes any quark and anti-quark up to the b quark. The latter was 
generated simultaneously with the initial-state QCD and final-state QED LLL subtractions activated 
for radiative processes which include a jet (g or g) in the final state. The initial- and final-state parton 
showers (PS) were fully applied down to Qo = 5.0 GeV in both event generations, where a backward- 
evolution PS was used for the initial-state PS. The forced photon-radiation PS supplemented with the 
FF radiations was applied to the fragmentation process with a constraint of E^ > 15 GeV, while the 
ordinary QCD PS was applied to the other processes. 

The generated events were passed to PYTHIA 6.425 [Ti] in order to simulate further small-Q^ 
phenomena down to the hadron level. The PYTHIA simulation was carried out with its default 
setting, except for parp(67) = 1.0 and parp(71) = 1.0, as usual [5]. The kinematical cuts in 
Eq. ([ij were applied to the hadron-level events. Although the PYTHIA simulation produces some 
soft photons, there is no ambiguity in the selection of the two candidate photons since the GR@PPA 
simulation produces only two photons. The event generation conditions were set very loose in order 
not to affect the final results. Note that the gg — >■ 77 and two-fragmentation processes are not included 
in the current simulation. 

The resultant Ai?(7-jet) distributions from the radiative processes are shown in Fig. |8l We obtain 
a smooth distribution plotted with filled circles by combining the distributions for the fragmentation 
process (solid histogram) and LLL-subtracted 95 — >■ 77 + g process (dashed histogram). In the frag- 
mentation process, Ai?(7-jet) has been derived from the parton- level qg — > 77-t-g events reconstructed 
as described in a previous section, while it has been derived directly from the generated parton-level 
event in the LLL-subtracted g5 — > 77 + g process. The combined result should be compared with the 
naive gg — >■ 77 + g ME-prediction shown with a dotted histogram, as in Fig. [2] We can see that they 
are in reasonable agreement in most of the regions, and that the divergence at A_R(7-jet) = in the 
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Figure 9: Ai?(7-jet) distribution of the radiative processes after the isolation cut is apphed. 
notation of the histograms and plot is the same as that in Fig. \8\ 



The 



naive ME-prediction is naturally regularized in the combined result. We can also see a substantial 
enhancement of the combined result with respect to the naive ME-prediction in a wide range. This 
enhancement can be attributed to an effect of multiple QCD radiations in the final state. Although the 
distribution is finite, the events are strongly concentrated in a small A_R(7-jet) region. The fragmen- 
tation photon is likely to be surrounded by hadrons when Ai?(7-jet) is small. Isolation requirements 
will suppress these events. On the other hand, isolation requirements will not be effective for other 
processes. The contribution from the (?9 — > 77 + 5 process is illustrated with a dot-dashed histogram 
in Fig. [51 This distribution does not have any concentration at small A_R(7-jet). 



We impose an isolation condition by using a cone Et which is defined at the hadron level as 



E' 



T,iso 



AB.<Ri^. 



Ej 



(18) 



where Et is the transverse component of the particle energy with respect to the beam direction. The 
sum is taken over all particles inside a given AR cone around the photon, excluding the photon itself 
and neutrinos. We take as i?iso = 0.4 and require that Et^iso < 15 GeV. The Ai?(7-jet) distribution 
after applying the cut is shown in Fig. [9l We can see that the isolation cut strongly suppresses the 
events in A_R(7-jet) < 0.4, as expected. The events in this region are dominantly produced by the 
fragmentation process. About 40% of the events from the fragmentation process were rejected by 
the isolation cut. On the other hand, the cut has a small impact on the LLL-subtracted radiative 
processes. Although, about 5% of the events were rejected, including the contribution from qq — )• 77, 
the reduction in the cross section is only 1% because the nearly equal number of positive- and negative- 
weight events exhibited at small Ai?(7-jet) in the LLL-subtracted gg — >■ 77 + q, as shown in Fig. [5J 

Figure [TU] shows the Et, iso distributions for the two photons in the fragmentation events. The 
distribution for the leading (highest- -Et) photon is shown with a dashed histogram, and that for 
the subleading photon with a solid histogram. The leading photon predominantly corresponds to the 
prompt photon produced by the qg -^ 75 interaction, and the subleading photons are mostly produced 
in the fragmentation of the quark. Parton showers produce low-energy partons in a wide area and 
the underlying-event simulation is activated in the PYTHIA simulation. These hadronic activities 
produce non-zero -Et.Iso even for the prompt photon by an accidental overlap. This contribution 
steeply decreases as Et, iso increases, as we can see in the dashed histogram in Fig. [TUl In contrast. 
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Figure 10: St.Iso distributions for tire selected two photons in the fragmentation events. The distribu- 
tion for the leading (highest-£'T) photon is shown with a dashed histogram and that for the subleading 
photon with a solid histogram. 

the distribution for the fragmentation photons shows a long tail extending to large i?T, iso- This tail 
corresponds to the hadronic activity produced in association with the photon, i.e., the remnant jet. 

Although small-Ait!(7-jet) events are strongly suppressed, the isolation cut does not produce a 
sharp edge in the Ai?(7-jet) distribution owing to a spread of the remnant jets. Figure [TT] shows 
the pt distribution of the remnant quark in the parton-level qg ^ ^^ + q events reconstructed for 
the fragmentation events. The distribution for those events inside the parton-level isolation cone, 
Ai?(7-jet) < 0.4, are plotted. The solid and dashed histograms show the distribution after and before 
the isolation cut, respectively. The events with PT(jet) > 15 GeV/c are likely to be rejected by the 
cut. However, the cut does not produce a sharp edge. On the other hand, a sharp edge should appear 
in this distribution and the distribution in Fig. ^ if we apply a similar cut at the parton level. This 
kind of smearing effect due to the final-state PS and hadronization cannot be evaluated without using 
simulations at the hadron level. 



6 Combined results 

We can obtain a consistent simulation result by combining the two simulations described in the 
previous section. Though we have not mentioned to the result from the qq — >■ 77 process previously, 
this process is simultaneously simulated with the LLL-subtracted radiative processes. The initial-state 
PS applied to the qq — >■ 77 events regularizes the initial-state LL components subtracted from the 
radiative processes. 

Figure [12] shows the pT distribution of the diphoton (77) system from the combined simulation. 
The isolation cut described in the previous section has been applied to the simulated events selected 
with the kinematical condition in Eq. ([T]). In addition to the combined result (solid histogram), the 
distributions of subprocesses are separately shown. We can see that the pT dependence is significantly 
different between the subprocesses. The integrated contribution from each subprocess is: 4.2 pb 
(23%) from qq -^ 77, 9.7 pb (54%) from the fragmentation, and 3.3 pb (18%) and 0.9 pb (5%) from 
the LLL-subtracted qg — >■ 77 -|- q and qq ^ ^^ + g processes, respectively. It is remarkable that the 
contribution from the lowest-order process, qq — >■ 77, is smaller than 1/4 of the sum. Furthermore, the 
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Figure 11: px distribution of the quark in the parton- level qg -^ 77 + q events reconstructed for the 
fragmentation events in Ai?(7-jet) < 0.4. The dashed histogram shows the distribution before the 
isolation cut and the solid histogram shows that after the cut, Sr.iso < 15 GeV. 
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Figure 12: px distribution of the diphoton (77) system in the simulated events satisfying the kinemat- 
ical condition in Eq. ([Ij and the isolation condition. The sum (solid histogram) is presented together 
with the distributions of the subprocesses: qq — ?> 77 (dashed histogram with dots), the LLL subtracted 
qg — > 77 + g (dashed histogram) and qq -^ jj + q (dot-dashed histogram), and the fragmentation 
process (dotted histogram). The result is compared with the predictions from RESBOS (filled circles) 
and DIPHOX (open circles). The DIPHOX prediction is presented at pT > 20 GeV/c only. 
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Pt(77) spectrum of qq — >■ 77 is apparently different from the summed spectrum. In the weak-boson 
production processes which we have previously studied, the production is totally dominated by the 
lowest-order process at small pt |13| . On the other hand, in the diphoton production, the lowest-order 
process is not dominant even at very small pt (^ 10 GeV/c) as we can see in Fig. [TH Studies with the 
lowest-order process alone are not sufficient even if lavge-pT contributions are intensively excluded. 
We can also see that simulations are not sufficient even if we include the fragmentation process in 
addition to the lowest-order process. The contribution from the LLL-subtracted qg — >■ 77 -t- q process, 
which supplements the fragmentation process in the hard-radiation region, is nearly comparable to 
the lowest-order process, and has a prill) spectrum remarkably different from the lowest-order and 
fragmentation processes. 

It must be noted that the composition described in the above is not physically meaningful. The 
above result is obtained when we separate the soft and hard radiations at pT of the qq — > 77 or 
qg — > 75 interaction. The composition changes if we adopt other definitions. Incidentally, it must 
also be noted that the fragmentation process is a part of the 95 — > 77 + q process. The fragmentation 
process has been introduced to regularize the final-state divergence in qg — )■ 77 -I- q. The sum of the 
fragmentation and LLL-subtracted qg — > 77 + q processes amounts to more than 70% of the total 
yield in our simulation. This fraction would not change dramatically even if we moved the separation 
criteria in a reasonable range. The domination of the 93 — >■ 77 + 9 process is a characteristic feature 
of the diphoton production in hadron collisions, at least, at the LHC with a typical Higgs-search 
condition. 

The simulated prill) spectrum is compared with the prediction from RESBOS [4] in Fig. [121 for 
which we have used CTEQ6M [15] for PDF and required the same kinematical and isolation conditions. 
RESBOS provides an NLO prediction in which initial-state soft QCD radiations are resummed. It is 
considered to be most reliable at present as far as the prill) distribution is concerned. Unfortunately, 
since the resummation can only predict inclusive properties, RESBOS cannot provide the hadron- 
level information. Thus, the selections are all applied at the parton level. The contribution from the 
99 ~^ 11 process is not included in the compared RESBOS prediction because this process is yet to be 
included in our simulation. The RESBOS prediction (filled circles) shown in Fig. [T^] is smaller than 
the result presented in the original paper [i] because the m^-y constraint in Eq. ([IJ is additionally 
required. We have confirmed that we can obtain a prediction close to the original result if the m-yy 
constraint is removed. As we can see in Fig. 1121 our simulation gives a substantially larger cross section 
than the RESBOS prediction although the overall tendency is consistent. The total cross section is 
15.6 pb from RESBOS and 18.1 pb from our simulation. 

The open circles in Fig. [H] show the prediction from DIPHOX [3] with the CTEQ6M PDF. 
DIPHOX is a full NLO calculation for diphoton production. The predictions are considered to be 
comparable with those from RESBOS. However, because soft QCD effects are not resummed, it 
cannot provide a continuous spectrum at small prill)- In Fig. [121 we have plotted the DIPHOX 
prediction in a relatively large prill) region (> 20 GeV/c) where we can obtain a continuous prill) 
prediction. We have required the kinematical condition in Eq. ([Ij and the same isolation condition. 
The plotted prediction is the sum of the results for the "direct" process and the one-fragmentation 
process separately evaluated. We have used the LO prediction for the one-fragmentation process, 
and subtracted the gg — >■ 77 contribution from the direct process. The obtained DIPHOX prediction 
is continuously smaller than our simulation and the RESBOS prediction in the displayed prill) 
range, although the difference becomes smaller as prill) becomes larger. This prill) behavior is 
reasonable because DIPHOX does not include initial-state radiation effects. The total cross section 
from this DIPHOX prediction is 13.7 pb. 

Figure [T3l compares the diphoton invariant-mass im,^^) distribution obtained from our simulation 
(solid histogram) with those from RESBOS (filled circles) and DIPHOX (open circles). The results are 
almost comparable to each other, although there are substantial overall shifts between them and small 
differences in the m^^ dependence. If we ignore the overall shifts, the DIPHOX prediction is closer 
to our simulation than RESBOS. The small-m-y^ fraction tends to be larger in our simulation and 
DIPHOX. Here, we have to remind that the isolation requirement in RESBOS is markedly different 
from that in DIPHOX and our simulation. Our simulation sums up the Er of particles inside the 
isolation cone, /S.R < 0.4, and rejects events if the summed Er, Er,iso, exceeds the threshold of 
15 GeV. In principle, DIPHOX does the same rejection at the parton level. On the other hand. 
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Figure 13: Invariant-mass {m-yy) distribution of the diphoton system. Together with the combined 
result (sohd histogram), the contribution from the fragmentation process (dotted histogram) and that 
from the other processes (dashed histogram) are separately presented. The result is compared with 
the predictions from RESBOS (filled circles) and DIPHOX (open circles). 

although the selection condition is literally identical, RESBOS rejects the whole LL contribution 
when Ai?(7-jet) < 0.4, in order to regularize the final-state divergence. Events with Et.iso < 15 GeV 
are allowed only for non-LL components. As we can see in Fig. |S1 the fragmentation process that 
simulates the LL component in our simulation dominate the cross section inside the isolation cone. 
The contribution from the fragmentation process and that from the other processes are separately 
shown with dotted and dashed histograms, respectively, in Fig.[T31 We can see that the fragmentation 
process shows a substantial enhancement at smaller m^^ with respect to the others, suggesting that 
the small-?7i^^ enhancement in DIPHOX and our simulation may be explained by the difference in 
the isolation requirement. The difference between our simulation and RESBOS, which we can see in 
Fig. 1121 seems to support this argument; the difference is large in the region where the fragmentation 
contribution is large. 

As a test, we have carried out a hybrid simulation in order to clarify the reason for the observed 
difference. In this simulation, we applied the kinematical cut and the isolation cut directly to the 
generated parton-level events in the LLL-subtracted radiative processes, without applying PS and 
the PYTHIA simulation. In the fragmentation process, only the final-state PS was applied to the 
generated events in order to reconstruct the parton-level qg — >■ 77 -I- q events. The kinematical cuts 
were applied to these parton-level events, and all events having Ai?(7-jet) < 0.4 were rejected in 
order to simulate the calculation in RESBOS. These are extensions of the simulations in the matching 
test. On the other hand, the full simulation down to the hadron level was applied to the lowest-order 
qq — >■ 77 process in order to make the combined result comparable with the resummed calculation in 
RESBOS. 

The combined Pt{ii) distribution from the hybrid simulation is presented with a solid histogram 
in Fig. [TH to be compared with the RESBOS prediction plotted with filled circles. We can see a 
complete agreement between them over almost the entire pt{ii) region, except for a narrow region 
at small prill)- The total cross section from the hybrid simulation is 16.5 pb, which is only 6% 
larger than the RESBOS prediction. The agreement is remarkable but not surprising because the 
radiative processes which dominate the total cross section are evaluated at the tree level even in 
NLO calculations, such as RESBOS and DIPHOX. The contribution of non-divergent soft/virtual 
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Figure 14: pT distribution of the diphoton (77) system obtained from the hybrid simulation. The sohd 
histogram shows the result from a RESBOS-style simulation in which all fragmentation events with 
Ai?(7-jet) < 0.4 are rejected. The dashed histogram has been obtained from a simulation in which 
the standard isolation cut, Etaso < 15 GeV, is applied at the parton level in all radiative processes. 
The results are compared with the RESBOS prediction shown with filled circles. 

corrections which are missing in our simulation must be small, since the contribution of the lowest- 
order process to which the corrections are to be applied is small. The difference at small px (77) may 
have been caused by the difference between the PS and the resummation. In any case, the agreement 
implies that the calculations are reasonably carried out both in RESBOS and our simulation. 

The dashed histogram in Fig. [T?] shows the result from another hybrid simulation, in which we 
allowed those events with i?T,iso < 15 GeV in the fragmentation process. This change resulted in a 
step- like structure in the pt{ii) spectrum. The reason for this is quite simple. The pt of the diphoton 
system is equal to the px of the remnant jet in the fragmentation process, and E't.Iso is equal to this px 
if the remnant jet is inside the isolation cone. Therefore, the allowance of non-zero Etaso produces an 
enhancement only at pt{ii) < 15 GeV. This change increased the total cross section to 18.1 pb, which 
is in good agreement with the full simulation result. This step-like enhancement can be considered 
as the source of the difference between our simulation and RESBOS in Fig. [T^l The enhancement 
must have been smeared by multiple QCD-radiation in PS and hadronization. Hence, the difference 
in Fig. [T^ must be reasonable and shows the effects that are not supported in RESBOS. 

We have carried out a study by using our custom-made PS in the above. We performed a similar 
study previously [20] by using the PYTHIA PS [M] for generating the fragmentation events. The 
matching properties were not fully satisfactory and the resultant diphoton-production cross section 
was smaller than the RESBOS prediction in that study. These were mainly caused by the fact that the 
photon yield from the PYTHIA PS was significantly smaller than that from our PS. The Ai?(7-jet) 
distributions obtained with the PYTHIA PS and our PS are compared in Fig. [151 We can see that 
the PYTHIA-PS predictions are almost always smaller than ours by about a factor of two. We do not 
understand the reason for this difference. In addition, though it may be too much detail, the drop 
in the first bin that we can see in the PYTHIA results is caused by the absence of small-Q^ photon 
radiations (Q < 1 GeV). Such small-Q^ radiations are simulated with the FF radiation in our PS. 
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Figure 15: Ai?(7-jet) distribution of the fragmentation process simulated with the PYTHIA PS and 
our PS in the standard setting. The distribution before and after the isolation cut is shown with a 
solid and dashed histograms, respectively, for the result with the PYTHIA PS. The corresponding 
distributions with our PS are plotted with filled circles and open circles, respectively. 

7 Conclusion 

We have developed an exclusive event generator for non-resonant QED diphoton (77) production at 
hadron collisions, consistently including additional one-jet production. The qg ^^ ^n + q process to be 
included as a radiative process has a final-state QED divergence together with an initial-state QCD 
divergence. We have developed a final-state subtraction method by extending the method developed 
for initial-state QCD divergences (the LLL subtraction). The subtraction works as well as we expected. 
The differential cross section becomes finite after the subtraction and converges to zero at the limit 
where the original cross section diverges. 

We have also developed a final-state parton shower (PS) for generating fragmentation events. The 
PS implements QED photon radiations as well as ordinary QCD radiations, and can enforce hard- 
photon radiation in order to improve the generation efficiency. Small-Q^ photon radiations which 
are not covered by the PS arc simulated by using a fragmentation function (FF). The divergent 
components in the final state which are subtracted from the radiative processes are restored with an 
appropriate regularization by applying this PS to qg — > 79 events. The subtraction and regularization 
of initial-state QCD divergences are carried out in the same way as applied to weak-boson productions. 
The generated fragmentation events show a complete matching with the LLL-subtracted qg — > 77 -I- q 
events. 

A consistent simulation sample has been obtained by combining the fragmentation events with 
the qq — >■ 77 and LLL-subtracted qq ^ JJ + g and qg — > 77 -I- q events simultaneously generated in 
GR@PPA. The event generation was tested for the LHC design condition, pp collisions at 14 TeV. A 
typical Higgs-boson search condition was required to the events simulated down to the hadron level. 
A typical isolation cut was also applied to the photons. We observed a reasonable suppression of the 
events in the colhnear region, together with a visible smearing in physical distributions due to the 
application of PS and hadronization. 

In our simulation, for which we have chosen all the energy scales to be equal to px of the non- 
radiative qq — > 77 or qg — > "fq events, the contribution from the lowest-order qq — > 77 process is smaller 
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than 1/4 of the total yield. The contribution from the LLL-subtracted qg ^ JJ + q process, which 
supplements the fragmentation process in hard-radiation regions, is comparable to it. A consistent 
inclusion of this process is necessary for reliable studies. If we combine the contributions from the 
fragmentation and LLL-subtracted qg -^ "fj + q processes, the total qg —> 774-9 contribution amounts 
to more than 70%. 

The combined simulation sample shows a behavior which reasonably agrees with the predictions 
from ncxt-to-lcading order (NLO) calculations by RESBOS and DIPHOX. Observed differences can 
be attributed to the different implementation of the isolation requirement, and higher-order QCD 
and hadronization effects which arc not supported in NLO calculations. It must be noted that the 
two-fragmentation and gg — > 77 and their higher orders are yet to be included in our simulation. 
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